---
title: "4 - PTSD Meta analysis"
author: "Alvaro Passi-Solar"
date: "`r Sys.Date()`"
output:
  word_document: default
  html_document: default
---


```{r, echo=FALSE, warning=FALSE, include=F}

library(weightr)
library(meta)
library(metasens)
library(metafor)
library(captioner)
library(tidyverse)

fig_nums<-captioner(prefix = "Figure 1.",  auto_space = FALSE)
tab_nums<-captioner(prefix = "Table 1.",  auto_space = FALSE)

settings.meta(digits = 3)
options(scipen=999)

FitFlextableToPage <- function(ft, pgwidth = 7){
  
  ft_out <- ft %>% flextable::autofit()
  
  ft_out <- flextable::width(ft_out, width = dim(ft_out)$widths*pgwidth /(flextable::flextable_dim(ft_out)$widths))
  return(ft_out)
}
```




```{r, echo=FALSE, warning=FALSE, include=F}
df0<-read.csv("df0_SR.csv")


```


```{r}
# birtcohortx<-"_birtcohort"
# birtcohortx<-""
if(birtcohortx=="_birtcohort"){
  df0_exc<-df0%>%
    subset(!is.na(vi))%>%
    subset(vi>0) %>%
    subset(grepl("Prevalence of comorbidities between mood and anxiety disorders",Title) | 
             grepl("Anxiety disorders in young people",Title) | grepl("Mental disorders and suicide risk in emerging",Title)| 
             grepl("Breastfeeding and mental health in adulthood",Title)| 
             grepl("Posttraumatic stress disorder: Prevalences",Title))
  
  
  df0<-df0%>%
    subset(!is.na(vi))%>%
    subset(vi>0) %>%
    subset(!grepl("Prevalence of comorbidities between mood and anxiety disorders",Title)) %>%
    subset(!grepl("Anxiety disorders in young people",Title)) %>%
    subset(!grepl("Mental disorders and suicide risk in emerging",Title)) %>%
    subset(!grepl("Breastfeeding and mental health in adulthood",Title)) %>%
    subset(!grepl("Posttraumatic stress disorder: Prevalences",Title))
}
```


```{r, include=T}
df0$Anyanx.noutcomes.<-df0$Anyanxnoutcomes
df0$anyx<-cut(df0$Anyanx.noutcomes.,c(2,5,6,7,8), labels = c("3-5","6","7","8"))
table(df0$anyx,df0$Anyanx.noutcomes, useNA = "always")


df0$MASIST_2<-df0$MASIST
df0$MASIST_2[df0$MASIST=="DSM-III-R"]<-"DSM"
df0$MASIST_2[df0$MASIST=="DSM-IV"]<-"DSM"
df0$MASIST_2[df0$MASIST=="DSM-V"]<-"DSM"
table(df0$MASIST_2,df0$MASIST)

# unique(df0$MASIST)
```





# Meta analysis


```{r ,  warning = FALSE, include=F}
# l1_o<-for (i in c("ANY1","ANY2","ANY4")){

fx1<-function(condx){
  # condx<-"PTSD1"
  print(condx)
  condx<-condx
  condx1<-substr(condx,nchar(condx),nchar(condx))
  condx0<-substr(condx,1,nchar(condx)-1)
  
  # subgx<-"Interview"
  fx2<-function(subgx){
    subgx<-subgx
    print(subgx)
    
    
    if (condx=="PTSD1"){
      df0_ANY1<-df0%>%
        subset(MAALL_j=="PTSD1" | MAALL_j=="99 (PTSD1)")
    }
    if (condx=="PTSD2"){
      df0_ANY1<-df0%>%
        subset(MAALL_j=="PTSD2" | MAALL_j=="99 (PTSD2)")
    }
    
    if (condx=="PTSD4"){
      df0_ANY1<-df0%>%
        subset(MAALL_j=="PTSD4" | MAALL_j=="99 (PTSD4)")
    }
    # }
    
    
    
    
    df0_ANY1$uno<-1
    df0_ANY1<-df0_ANY1 %>%
      group_by(Estudio, label)%>%
      mutate(ticker=cumsum(uno))
    
    df0_ANY1<-df0_ANY1 %>%
      subset(ticker==1 | (ticker>1 & grepl("Medina-Mora",ID)))
    
    df0_ANY1$Interview_wo_dis<-df0_ANY1$Interview
    df0_ANY1$Interview_wo_dis[df0_ANY1$Interview=="DIS"]<-NA
    
    sg_f<-function(subgroupx) {
      
      df0_ANY1$subgroup<-df0_ANY1[[subgx]]
      
      df0_ANY1<-df0_ANY1 %>%
        subset(!is.na(subgroup))
      
      df0_ANY1<-df0_ANY1 %>%
        group_by(Estudio, subgroup)%>%
        mutate(sum_=sum(uno))
      
      df0_ANY1<-df0_ANY1 %>%
        select(sum_,everything()) %>%
        subset(sum_==1 | (sum_>1 & (!grepl("99",MAALL_j) )))
      
      
      
      rma_f<-function(subx){
        
        print(subx)
        d<-rma(yi=yi,vi=vi,
               data=subset(df0_ANY1,subgroup==subx),
               digits=5,
        )
        
        p<-predict(d,transf=transf.ilogit,
                   # targ=list(ni=df0_ANY1$total)
        )
        # p$coef<-coef(d)
        # p$se<-coef(d)
        p$subgroup<-subx
        p<-as.data.frame(p)
        p$estimate<-d$b[1]
        p$stderror<-d$se
        p$tau2s<-round(d$tau2,2)
        # p$tau2<-res$tau2
        p$Q<-round(d$QE,2)
        p$df<-d$k - d$p
        p$p<-metafor:::.pval(d$QEp,showeq=TRUE)
        p$I2<-round(d$I2,1)
        # p$I2<-res$I2
        # p$text<-mlabfun("",d)
        p
      }
      
      
      df0_ANY1$subgroup<-df0_ANY1[[subgroupx]]
      if(length(unique(df0_ANY1$subgroup))!=1){
        pred_m1_f<-function(){
          m1<-rma(yi=yi,vi=vi,
                  data=df0_ANY1,
                  method="REML",
                  slab=ID,
                  random = ~ 1 | ID,
                  digits=5)
          pred_m1<-predict.rma(m1, transf=transf.ilogit)
          
          pred_m1<-pred_m1 %>%
            as.data.frame()
          pred_m1$estimate<-m1$b[1]
          pred_m1$stderror<-m1$se
          pred_m1$tau2s<-round(m1$tau2,2)
          # p$tau2<-res$tau2
          pred_m1$Q<-round(m1$QE,2)
          pred_m1$df<-m1$k - m1$p
          pred_m1$p<-metafor:::.pval(m1$QEp,showeq=TRUE)
          pred_m1$I2<-round(m1$I2,1)
          pred_m1$text<-paste0("Heterogeneity: ",
                               " I²=", pred_m1$I2,"%",
                               ", T²=",pred_m1$tau2s,
                               ", Q=",pred_m1$Q,
                               ", df=",pred_m1$df,
                               " (",
                               "p",pred_m1$p,
                               ")")
          list(m1,pred_m1)
        }
        pred_m1<-pred_m1_f()[[2]]
        pred_m1<-pred_m1%>%
          as.data.frame()
        pred_m1$subgroup<-"Total"
        pred_m1$ID_AP<-"RE Total"
        total_tex_pi<-paste0(
          # "",
          #                format(round(adf$pred*100, digits=2), nsmall = 2),
          "PI: ",
          format(round(pred_m1$pi.lb*100, digits=2), nsmall = 2),
          "-",
          format(round(pred_m1$pi.ub*100, digits=2), nsmall = 2),"")
        
        pred_m1b<-pred_m1
        pred_m1b$subgroup<-"Total"
        pred_m1b$ID_AP<-""
        pred_m1b$pred<-NA
        pred_m1b$ci.lb<-NA
        pred_m1b$ci.ub<-NA
        
        pred_m1$text<-""
        pred_m1<-bind_rows(pred_m1,pred_m1b)
        
        adf_fl<-lapply(unique(df0_ANY1$subgroup),
                       rma_f)
        
        adf<-do.call(bind_rows,adf_fl)
        adf<-adf%>%
          select(subgroup,everything())
        # as.character(adf$text)
        adf$text<-paste0("Heterogeneity: ",
                         " I²=", adf$I2,"%",
                         ", T²=",adf$tau2s,
                         ", Q=",adf$Q,#X²
                         # \n
                         ", df=",adf$df,
                         " (",
                         "p",adf$p,
                         ")")
        
        adf$tex_pi<-paste0(
          # "",
          #                format(round(adf$pred*100, digits=2), nsmall = 2),
          "PI: ",
          format(round(adf$pi.lb*100, digits=2), nsmall = 2),
          "-",
          format(round(adf$pi.ub*100, digits=2), nsmall = 2),"")
        
        res<-rma.mv(yi, vi, 
                    method="REML",
                    mods = ~ subgroup, 
                    random = ~ subgroup | ID, 
                    slab=ID,
                    # struct="ID", 
                    data=df0_ANY1, 
                    digits=5)
        
        
        adf$overall_QM<-res$QM
        adf$overall_df<-res$p-1
        adf$overall_QMp<-res$QMp
        
        pd<-df0_ANY1 %>%
          tibble()%>%
          arrange(Country, Fieldwordmidpoint)
        
        # adf$Country<-adf$subgroup
        adf0<-adf %>%
          select(subgroup,text,tex_pi,overall_QM,overall_df,overall_QMp)
        adf0$ID_AP<-""
        
        dd<-adf %>%
          select(subgroup,pred,ci.lb,ci.ub)
        dd$ID_AP<-"RE subgroup"
        
        
        dd_adf0<-bind_rows(dd,adf0)
        
        pd<-bind_rows(pd,dd_adf0)
        
        
        
        
        
        pd<-bind_rows(pd,pred_m1)
        
        pd<-pd %>%
          arrange(ID_AP)
        
        
        pd$label_ct<-paste0(
          round(pd$cases,0),
          " (",
          round(pd$total,0),
          ")")
        
        pd$label_ct[pd$ID_AP==""]<-""
        pd$label_ct[pd$ID_AP=="Total"]<-""
        pd$label_ct[pd$ID_AP=="RE subgroup"]<-""
        pd$label_ct[pd$ID_AP=="RE Total"]<-""
        
        pd$label<-paste0(format(round(pd$P*100, digits=2), nsmall = 2),
                         " (",
                         format(round(pd$inf_p*100, digits=2), nsmall = 2),
                         "-",
                         format(round(pd$sup_p*100, digits=2), nsmall = 2),
                         ")")
        
        pd$label[is.na(pd$P)]<-paste0(format(round(pd$pred[is.na(pd$P)]*100, digits=2), nsmall = 2),
                                      " (",
                                      format(round(pd$ci.lb[is.na(pd$P)]*100, digits=2), nsmall = 2),
                                      "-",
                                      format(round(pd$ci.ub[is.na(pd$P)]*100, digits=2), nsmall = 2),
                                      ")")
        
        pd$label[pd$ID_AP==""]<-""
        
        
        
        pd$size<-1/ sqrt(pd$vi)
        pd$size <- (pd$size/20) / max(pd$size, na.rm = T)
        pd$ID_AP[!is.na(pd$text)]<-pd$text[!is.na(pd$text)]
        
        pd$text_caption[pd$ID_AP!=""]<-unique(pd$text_caption[pd$ID_AP==""])
        pd$P[!is.na(pd$pred)]<-pd$pred[!is.na(pd$pred)]
        pd$inf_p[!is.na(pd$ci.ub)]<-pd$ci.ub[!is.na(pd$ci.ub)]
        pd$sup_p[!is.na(pd$ci.lb)]<-pd$ci.ub[!is.na(pd$ci.lb)]
        pd$pred_t_l<-pred_m1$ci.lb[1]
        pd$pred_t<-pred_m1$pred[1]
        pd$pred_t_u<-pred_m1$ci.ub[1]
        pd$nchar<-nchar(pd$label_ct)
        maxlabelct<-max(pd$nchar)
        pd$nchar_max<-maxlabelct-pd$nchar
        
        
        for(i in 1:length(pd$label_ct)){
          pd$label_ct_t[i]<-paste0(pd$label_ct[i],paste0(rep(" ", times=1+pd$nchar_max[i]), collapse = ""))
        }
        pd$nchar<-nchar(pd$label)
        maxlabelct<-max(pd$nchar)
        pd$nchar_max<-maxlabelct-pd$nchar
        
        pd$label[!is.na(pd$tex_pi)]<-pd$tex_pi[!is.na(pd$tex_pi)]
        
        pd$label[pd$label=="" & pd$subgroup=="Total" & pd$ID_AP!="Total"]<-total_tex_pi
        
        
        for(i in 1:length(pd$label_ct)){
          pd$label_2[i]<-paste0(" ",pd$label[i],paste0(rep(" ", times=pd$nchar_max[i]), collapse = ""))
        }
        # nchar(pd$label_ct_t)
        
        pd$label_concat<-paste0(pd$label_ct_t,pd$label_2)
        # pd$label_concat[pd$ID_AP=="RE subgroup"]<-paste0("         ",pd$label_concat[pd$ID_AP=="RE subgroup"])
        pd$label_concat<-gsub("NA","", pd$label_concat)
        # nchar(pd$label_concat)
        pd_sub<-pd %>%
          subset(!is.na(subgroup))%>%
          subset(subgroup!="Total")%>%
          select(subgroup)%>%
          unique()
        pd_sub$order<-0
        pd_sub$Fieldwordmidpoint<-0
        pd_sub$ID_AP<-pd_sub$subgroup
        pd_sub$bold<-1
        
        pd<-bind_rows(pd,pd_sub)
        pd<-pd %>%
          select(order, everything())
        
        
        pd$ID_AP[pd$subgroup=="Total" & pd$ID_AP==""]<-"Total"
        pd$order[pd$ID_AP=="Total"]<-7777
        pd$Fieldwordmidpoint[pd$ID_AP=="Total" ]<-7777
        
        pd$order[grepl("RE s",pd$ID_AP)]<-8888
        pd$order[grepl("Hete",pd$ID_AP)]<-9999
        
        pd$Fieldwordmidpoint[grepl("RE s",pd$ID_AP)]<-8888
        pd$Fieldwordmidpoint[grepl("Hete",pd$ID_AP)]<-9999
        
        pd<-pd %>%
          arrange(subgroup,Fieldwordmidpoint,order)
        # pd$Fieldwordmidpoint
        
        pd$order_final<-1:nrow(pd)
        
        pd$nchar<-nchar(pd$ID_AP)
        maxlabelct<-max(pd$nchar)
        pd$nchar_max<-maxlabelct-pd$nchar
        
        for(i in 1:length(pd$ID_AP)){
          pd$ID_AP[i]<-paste(pd$ID_AP[i],paste0(rep(" ",times=pd$nchar_max[i]),collapse = ""))
        }
        # nchar(as.character(pd$ID_AP))
        # pd$ID_AP<-paste(pd$ID_AP,paste0(rep(" ",max(nchar(pd$ID_AP))-nchar(pd$ID_AP)),collapse = ""))
        pd$label_concat[is.na(pd$label_concat)][1]<-"n (N)        %  (95% CI)"
        pd$label_concat[is.na(pd$label_concat)]<-""
        
        pd$ID_AP<-paste0(pd$ID_AP,pd$label_concat)
        pd$ID_AP<-gsub("( ","(",pd$ID_AP,fixed=T)
        pd$ID_AP<-gsub("- ","-",pd$ID_AP,fixed=T)
        pd$nchar<-nchar(pd$ID_AP)
        maxlabelct<-max(pd$nchar)
        pd$nchar_max<-maxlabelct-pd$nchar
        
        for(i in 1:length(pd$ID_AP)){
          pd$ID_AP[i]<-paste(pd$ID_AP[i],paste0(rep(" ",times=pd$nchar_max[i]),collapse = ""))
        }
        
        
        pd<-pd %>%
          arrange(order_final)%>%
          mutate(ID_AP=factor(ID_AP, levels=unique(ID_AP)) )%>%
          select(order_final,everything())
        
        
        pd<-pd %>%
          select(order_final, everything())%>%
          arrange(order_final)
        pd$y_disc<-factor(pd$order_final)
        pd$y_disc<-fct_rev(pd$y_disc)
        pd<-pd %>%
          select(y_disc, everything())
        
        pd$label_concat<-gsub("NA","",pd$label_concat)
        pd$face<-"plain"
        pd$face[pd$bold==1 & !is.na(pd$bold)]<-"bold"
        pd$face[grepl("Total",pd$ID_AP)]<-"bold"
        
        pd$ID_AP<-fct_rev(pd$ID_AP)
        
        # ))
        pd$text_caption<-paste0("Test of Moderators QM (df=",
                                round(res$p-1,0),
                                ") = ",
                                round(res$QM,4),
                                ", p-val = ",
                                format(round(res$QMp,4), nsmall=4))
        pd<-pd%>%
          select(face, everything())
        
        # windows()
        pd$ID_AP<-sub("         ","",pd$ID_AP,fixed=T)
        
        # showtext::showtext_auto()
        
        p1<-ggplot(pd) +
          geom_vline(aes(xintercept = pred_t_l), col="#3A86FF")+
          geom_vline(aes(xintercept = pred_t), col="#3A86FF")+
          geom_vline(aes(xintercept = pred_t_u), col="#3A86FF")+
          
          
          geom_point(aes(x=P, y=y_disc, size=size), shape=15, col="#FFBE0B")+
          geom_errorbarh(aes(xmax = sup_p, xmin = inf_p, 
                             y=y_disc),  height = .2, col="#FFBE0B")+
          scale_y_discrete(
            breaks= pd$y_disc,
            labels = pd$ID_AP
          )+
          
          geom_point(aes(x=pred, y=y_disc), shape=1, col="#FF006E")+
          geom_errorbarh(aes(xmax = ci.ub, xmin = ci.lb, y=y_disc),  height = .2, col="#FF006E")+
          scale_x_continuous(labels = percent)+
          # coord_cartesian(xlim=c(0,max(pd$sup_p, na.rm = T)*1.4))+
          # theme_f()+
          theme(
            text=element_text(family="mono"),
            axis.ticks.y = element_blank(),
            legend.position ="none",
            # aspect.ratio = 1,
            panel.background = element_rect(fill="white"),
            strip.background = element_blank(),
            strip.placement = "outside",
            strip.text = element_blank(),
            axis.text.y = element_text(size=10, colour = "black", angle = 0, hjust = 0, vjust = 0.5,
                                       face = pd$face),
            axis.text.x = element_text(size=10, colour = "black", angle = 0, hjust = 0.5, vjust = 0.5),
            axis.line.x=element_line(color="black"),
            plot.caption.position = "plot",
            plot.caption = element_text(hjust = 0)
          )+
          labs(
            title=condx,
            caption=unique(pd$text_caption),
            y="",
            x="")
        # /home/passi/Dropbox/
        ggsave(file=paste0("Render/4 - ",condx,"_forest_plot_",subgx,birtcohortx,".pdf"), plot=p1, width=15, height=10, dpi=300)
        ggsave(file=paste0("Render/png/4 - ",condx,"_forest_plot_",subgx,birtcohortx,".png"), plot=p1, width=15, height=10, dpi=600)
        list(p1,pd)
      }
      # table(pd$face)
    }
    sg_f(subgx)
    
  }
  
  lapply(c("MASIST_2","Interview","Interview_wo_dis"),
         fx2)
}
```


```{r ,  warning = FALSE, include=F}
l1o<-lapply(c("PTSD1", "PTSD2","PTSD4"),
            fx1)




```


```{r ,  warning = FALSE, include=F}
l1o
```




